rm(list=ls(all=TRUE))
set.seed(610)
library(foreign)
library(xlsx)
options(stringsAsFactors = FALSE)
Sys.setlocale(,"CHT")
deathorigin=read.csv("C:/Users/weary/Dropbox/Duke Prepare/11 the Final!/white terror studies project/repression list/final clean.csv")
head(deathorigin)
table(deathorigin$number)
deathorigin$taiwan[is.na(deathorigin$taiwan)]=0
table(deathorigin$taiwan)
table(deathorigin$area_origin3_clean,deathorigin$taiwan)
final=read.csv("C:/Users/weary/Dropbox/Duke Prepare/11 the Final!/white terror studies project/repression list/election for combine_after2000.csv")


table(deathorigin$area_origin3_clean,deathorigin$taiwan)



final$killbefore=NA
final$killbeforeTW=NA
for(i in 1:nrow(final)){
  test1=deathorigin[deathorigin$area_origin3_clean==final$county_name[i],]
  final$killbefore[i]=sum(test1$Clear_year3<=final$year[i]&test1$Clear_year3>final$year[i]-4,na.rm=T)
  final$killbeforeTW[i]=sum(test1$Clear_year3<=final$year[i]&test1$Clear_year3>final$year[i]-4&test1$taiwan==1,na.rm=T)
}

final$killbeforeS=NA
final$killbeforeTWS=NA
for(i in 1:nrow(final)){
  test1=deathorigin[deathorigin$area_origin3_clean==final$county_name[i],]
  final$killbeforeS[i]=sum(test1$Clear_year3<=final$year[i],na.rm=T)
  final$killbeforeTWS[i]=sum(test1$Clear_year3<=final$year[i]&test1$taiwan==1,na.rm=T)
}

final$killmd= final$killbefore-final$killbeforeTW
head(final)
final$group3=NA
final$group3[final$KMT_per>=0.5]="KMT"
final$group3[final$KMT_per<=0.5]="DunWai"


final$halfKMT=NA
final$halfKMT[final$KMT_last<=0.5]=0
final$halfKMT[final$KMT_last>=0.5]=1

final$perkill=final$killbeforeTWS/final$number_total_vote


#====Figure 1, Left Panel, Number of Victim per year======#

ydata=NULL
ydata$year=seq(1947,2016,1)
ydata=data.frame(ydata)
head(ydata)
ydata$allkill=NA
ydata$KMTrate=NA

for(i in 1:nrow(ydata)){
  ydata$allkill[i]=sum(deathorigin$Clear_year3==ydata$year[i],na.rm=T)  
  ydata$KMTrate[i]=sum(final$number_KMT_vote[final$year==ydata$year[i]],na.rm=T)/sum(final$number_total_vote[final$year==ydata$year[i]],na.rm=T)  
}

ggplot(data=ydata,aes(x=year,y=allkill))+
  geom_line()+geom_point(size=1)+theme_bw()+
  scale_x_continuous(breaks=seq(1940,2020,10))+ 
  ylab("Number of Victims")+ 
  geom_vline(xintercept=1987,lty="dashed")+ 
  annotate("text", x = 1987, y = 1000, label = "1987 \n Lift of Martial Law")

ggsave("Figure 1 L.png",width=5,height=3,units="in",dpi=300)

#===============================================================#
#Figure 3 and Regression


table(deathorigin$area_origin3_clean,deathorigin$taiwan)


combine2=NULL
combine2$chinese_name2 <- c("������", "�s����", "�s�˥�", "������","������", 
                            "�x�_��", "�x����", "�x����", "�x�n��","�x�n��", 
                            "�x�_��","��鿤", "���ƿ�", "�Ÿq��", "�Ÿq��",
                            "�s�˿�", "�Ὤ��", "�򶩥�", "�]�߿�",
                            "�n�뿤", "���", "�̪F��", "�x�F��",
                            "�y����", "���L��")

combine2$kill1987tw <- c(0,0,20,237,286,
                         518,161,406,187,456,
                         398,448,203,60,259,
                         254,33,61,337,
                         122,45,198,18,
                         135,208)


combine2$kill1987 <- c(241,130,39,826,371,
                       653,238,484,267,487,
                       1020,561,259,118,295,
                       102,155,213,356,
                       194,168,266,105,
                       168,243)



combine2$population1990 <- c(42754,5585,324426,1386723,1119263,
                             3048034,761802,1258157,683251,1026983,
                             2719659,1355175,1245288,257597,552277,
                             374492,352233,352919,547609,
                             536479,95932,893282,256803,
                             450943,753639)


combine2$population1966 <- c(56842,15868,267439,632662,732773,
                             1123354,380505,704290,416009,884817,
                             1331008,609979,991538,257597,552277,
                             267439,307220,287156,487317,
                             475315,112852,760101,267336,
                             384420,763423)


combine2$dens=combine2$kill1987/combine2$population1990
combine2$denstw=combine2$kill1987tw/combine2$population1990
combine2$dens2=combine2$kill1987/combine2$population1966
combine2$denstw2=combine2$kill1987tw/combine2$population1966

combine2=data.frame(combine2)

final$perkill3=NA
final$perkill3tw=NA

for(i in 1:nrow(final)){
  final$perkill3[i]= combine2$dens[combine2$chinese_name2==final$county_name[i]] 
  final$perkill3tw[i]= combine2$denstw[combine2$chinese_name2==final$county_name[i]] 
  final$perkill32[i]= combine2$dens2[combine2$chinese_name2==final$county_name[i]] 
  final$perkill3tw2[i]= combine2$denstw2[combine2$chinese_name2==final$county_name[i]] 
}

final$group3=NA
final$group3[final$KMT_per>=0.5]="KMT>50%"
final$group3[final$KMT_per<=0.5]="<50%"

ggplot(data=final[final$year>=1987&final$year<2014&final$perkill3<0.001,],aes(x=perkill3,y=KMT_per,colour=group3,fill=group3,size=number_total_vote,weight=number_total_vote))+
  geom_smooth(method = "lm", se = TRUE, show.legend = F)+geom_point()+theme_bw()+
  scale_colour_manual(values = c("green","blue"))+
  scale_fill_manual(values = c("green","blue"), guide= F)+
  labs(title="", x ="(N of Victims/Population1990)", 
       y= "KMT voteshare in county elections",
       size="N of Voters",colour="")
#final$perkill3<0.001 means to remove Kinmen, Matsu, and Penghu (p.16)
ggsave("Figure 3.png",width=8,height=6,units="in",dpi=300)

#=================================================================#
#Regression Table 1 (p.20)


final$KMTpres=0
final$KMTpres[final$year<2000]=1
final$KMTpres[final$year>=2008&final$year<=2015]=1


final$halfKMT=NA
final$halfKMT[final$KMT_per<=0.5]=0
final$halfKMT[final$KMT_per>=0.5]=1

final$lastKMT=NA
final$lastKMT[final$KMT_last<=0.5]=0
final$lastKMT[final$KMT_last>=0.5]=1


final$stableKMT=NA
final$stableKMT[final$KMT_last<=0.5&final$KMT_last2<=0.5]=0
final$stableKMT[final$KMT_last>=0.5&final$KMT_last2>=0.5]=1


model1=lm(KMT_per~perkill3*halfKMT,data=final[final$year>=1987&final$perkill3<0.001&final$year<2014,],weight=number_total_vote)
summary(model1)

model11=lm(KMT_per~perkill3*halfKMT+KMT_last+KMTpres,data=final[final$year>=1987&final$perkill3<0.001&final$year<2014,],weight=number_total_vote)
summary(model11)

model2=lm(KMT_per~perkill3*lastKMT,data=final[final$year>=1987&final$perkill3<0.001&final$year<2014,],weight=number_total_vote)
summary(model2)

final$perkill3.sd=final$perkill3-mean(final$perkill3,na.rm=T)

model22=lm(KMT_per~perkill3*lastKMT+KMTpres,data=final[final$year>=1987&final$perkill3<0.001&final$year<2014,],weight=number_total_vote)
summary(model22)

vif(model22)

model3=lm(KMT_per~perkill3*stableKMT,data=final[final$year>=1987&final$perkill3<0.001&final$year<2014,],weight=number_total_vote)
summary(model3)

model33=lm(KMT_per~perkill3*stableKMT+KMT_last+KMTpres,data=final[final$year>=1987&final$perkill3<0.001&final$year<2014,],weight=number_total_vote)
summary(model33)

library(stargazer)
stargazer(model1,model11,model2,model22,model3,model33,type="text",omit.stat=c("LL","ser","f"))

#======#Figure 1 R p.15====================================================#


table(deathorigin$area_origin3_clean)
require("knitr")
require("kableExtra")
require("maptools")
require("plyr")
require("ggplot2")
library("rgdal")
taiwan_shp <- readShapeSpatial("C://Users//weary//Dropbox//Duke Prepare//11 the Final!//white terror studies project//repression list//TWN_adm2.shp")
print(as.character(taiwan_shp$NAME_2))
taiwan_map <- fortify(taiwan_shp)
ggplot(taiwan_map, aes(x = long, y = lat, group=group)) +
  geom_path() + 
  coord_map() +xlim(118,122)+ylim(21.5,26.2)


mydata <- data.frame(NAME_1=taiwan_shp$NAME_2,
                     id=taiwan_shp$ID_2)

mydata$population1990 <- c(42754,5585,324426,2505986,
                           3048034,2019959,1710234,2719659,
                           1355175,1245288,257597,552277,
                           374492,352233,352919,547609,
                           536479,95932,893282,256803,
                           450943,753639)

mydata$kill1987 <- c(241,130,39,1197,
                     653,722,757,1020,
                     561,259,118,295,
                     102,128,213,356,
                     194,168,266,105,
                     168,243)

taiwan_map$id <- as.character(as.integer(taiwan_map$id)+1)
final.plot<-merge(taiwan_map,mydata,by="id",all.x=T)
head(final.plot)
library(RColorBrewer) 
library(mapproj)
final.plot$dens=final.plot$kill1987/final.plot$population1990

ggplot() +
  geom_polygon(data = final.plot[final.plot$dens<=0.001,],  #<0.001 for removing the three islands
               aes(x = long, y = lat, group = group, 
                   fill = kill1987), 
               color = "black", size = 0.25) + theme_bw()+
  coord_map()+
  scale_fill_gradientn(colours = brewer.pal(9,"Reds"), name = "N of Victims")+
  xlim(119.5,122)+ylim(21.9,25.3)

ggsave("Figure 1 R.png",height=4,width=4,units="in",dpi=300)

#======#Figure 2  (p.16)===========================#
ggplot() +
  geom_polygon(data = final.plot[final.plot$dens<=0.001,], 
               aes(x = long, y = lat, group = group, 
                   fill = kill1987/population1990), 
               color = "black", size = 0.25) + theme_bw()+
  coord_map()+
  scale_fill_gradientn(colours = brewer.pal(9,"Reds"), name = "(N of Victims)/\n(Population 1990)")+
  xlim(119.5,122)+ylim(21.9,25.3)+
  labs(title="", x ="", y = "")

ggsave("Figure 2.png",height=6,width=6,units="in",dpi=300)
